Rare-gas solids under pressure: A path-integral Monte Carlo simulation 
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Rare-gas solids (Ne, Ar, Kr, and Xe) under hydrostatic pressure up to 30 kbar have been stud- 
ied by path-integral Monte Carlo simulations in the isothermal-isobaric ensemble. Results of these 
simulations have been compared with available experimental data and with those obtained from a 
quasiharmonic approximation (QHA). This comparison allows us to quantify the overall anharmonic- 
ity of the lattice vibrations and its influence on several structural and thermodynamic properties of 
rare-gas solids. The vibrational energy increases with pressure, but this increase is slower than that 
of the elastic energy, which dominates at high pressures. In the PIMC simulations, the vibrational 
kinetic energy is found to be larger than the corresponding potential energy, and the relative differ- 
ence between both energies decreases as the applied pressure is raised. The accuracy of the QHA 
increases for rising pressure. 

PACS numbers: 67.80.-s, 62.50.+p, 65.40.De, 05.10.Ln 
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I. INTRODUCTION 

The importance of anharmonic effects in solids has 
been recognized long ago, as they are responsible for well- 
known phenomena such as thermal expansion, pressure 
dependence of the compressibility, phonon couplings, as 
well as isotope dependence of structural properties and 
melting temperature.^'^ This kind of effects have been 
studied both theoretically and experimentally for rare- 
gas solids since many years;^ because these are sim- 
ple systems allowing fruitful comparisons between theory 
and experiment. The interatomic forces are weak, short 
range, and rather well understood, so that critical tests 
of appropriate theories by their ability to predict prop- 
erties of actual rare-gas crystals are relatively simple. In 
particular, their thermodynamic properties are interest- 
ing due to the large anharmonic contributions to their 
lattice dynamics. 

From a theoretical point of view, anharmonic effects 
in solids have been traditionally studied by using ap- 
proaches such as the so-called quasiharmonic approxi- 
mation (QHA).^'^ In this approach, frequencies of vibra- 
tional modes are assumed to change with crystal vol- 
ume, and for given volume and temperature, the solid 
is supposed to be harmonic. ^'^ However, the QHA does 
not deal with phonon interaction effects, which can be 
treated by perturbation theorji when anharmonicity is 
not large, or by different self-consistent phonon theories 
for larger anharmonicitiesA^iiSiii A different theoretical 
procedure is the Feynman path integral methodf-^ which 
is well-suited to study thermodynamic properties of solids 
at temperatures lower than the Debye temperature 6_d, 
where the quantum nature of the atomic nuclei is rel- 
evant. The combination of path integrals with Monte 
Carlo (MC) sampling enables us to carry out quantita- 
tive and nonperturbative studies of anharmonic effects in 
solids. The path-integral Monte Carlo (PIMC) technique 
has been applied earlier to study several properties of 
rare-gas solids li^ii^ii^iiSiiLiSiiS In particular, it has pre- 



dicted kinetic-energy values in good agreement with ex- 
perimental dataiSSiSi An effective-potential Monte Carlo 



theoi 



has been also applied to study thermal and 



elastic properties of solid neon. 

Anharmonic effects increase appreciably with tem- 
perature. This is now well-known and has been ex- 
plained quantitatively for rare-gas solids. "'^'^^ In recent 
years, the effect of pressure on these solids has attracted 
much attention from both experimentalists^iSSiSSiSl and 
theorists iS'i^'SSiSS The influence of pressure on the an- 
harmonicity of lattice vibrations is, however, not well 
understood. It has been recently suggested that pres- 
sure causes a decrease in this anharmonicity, '^''■^'^ in line 
with earlier observations that the accuracy of the QHA 
increases as pressure is raisedi^ It has been also ar- 
gued that at high pressures, thermodynamic properties of 
solids can be well described by classical calculations, i.e., 
dealing with the atoms as classical oscillators in a given 
potentialiS This seems to be at first sight contradictory 
with the fact that pressure induces a larger zero-point 
vibrational energy of the solid. These questions are in- 
deed related with the ratio of the vibrational energy to 
the whole internal energy on one side, and with the size 
of the "intrinsic" anharmonicity (further than the QHA) 
of the lattice vibrations, on the other side. 

In this paper, we study structural and thermodynamic 
properties of rare-gas solids under pressure. This allows 
us to study properties of these solids along well-defined 
isotherms, and to analyze changes in anharmonic effects 
due to the repulsive (for compression) and attractive (for 
dilation, i.e., negative pressure) parts of the interatomic 
potential. The interatomic interaction is described by 
a Lennard- Jones potential. Results of the PIMC simu- 
lations are compared with those yielded by a quasihar- 
monic approximation with the same interatomic poten- 
tial. This approach will help us to quantify the influence 
of the "intrinsic" anharmonicity on the considered prop- 
erties. 

The paper is organized as follows. In Sec. II, the com- 
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TABLE I: Parameters a and e of the Lennard- Jones potential employed in this work and average isotopic mass (Af) for rare 
gases. Calculated zero-temperature properties of rare-gas solids at zero pressure are also given: lattice parameter a, zero-point 
vibrational energy E^ih, and elastic energy Ed per atom, as derived from PIMC simulations. 

Element a (A) e (meV) (M) (amu) a (A) E^a, (meV) Ed (meV) 



Ne 


2.782 


3.084 


20.18 


4.4631 


6.33 


1.20 


Ar 


3.404 


10.32 


39.95 


5.3115 


7.99 


0.43 


Kr 


3.638 


14.17 


83.80 


5.6458 


6.27 


0.18 


Xe 


3.961 


19.91 


131.30 


6.1316 


5.54 


0.10 



putational method is described. In Sec. Ill, we present 
results for energy, heat capacity, lattice parameter, and 
bulk modulus. Finally, Sec. IV includes a discussion of 
the results and the conclusions. 



II. METHOD 
A. Path-integral Monte Carlo 

Rare-gas atoms were treated as quantum particles in- 
teracting through a Lennard- Jones potential: V{r) = 
4e[(cr/r)^^ — ((j/r)^], with parameters e and a given in Ta- 
bic I, which were employed in earlier simulations of this 
kind of crystalspiSiiS In this Table we also give the average 
atomic mass of rare gases used in the calculations, as well 
as low-temperature properties of the studied crystals at 
zero pressure, as derived from PIMC simulations (see be- 
low) . Lennard- Jones-type potentials have been employed 
in recent years to model the atomic interaction in rare- 
gas solids iiSiiL22i2^ Although more sophisticated inter- 
action potentials have been developed, they do not seem 
to be significantly superior to Lennard- Jones potentials 
in accounting for the experimental data.^^'^^ This is not 
the case when one considers rare-gas solids under high 
pressure, where three-body potentials are necessary (see 
below) . For this reason, our calculations are restricted to 
pressures not higher than 30 kbar. 

Equilibrium properties of rare-gas solids have been cal- 
culated by PIMC simulations in the isothcrmal-isobaric 
ensemble (NPT). Simulations have been performed on 
5x5x5 cubic supercells of the face-centered-cubic unit 
cell, including 500 rare-gas atoms, and assuming peri- 
odic boundary conditions. To check the convergence of 
our results with system size, some MC simulations were 
carried out for other supercell sizes, including 7x7x7 
supercells. We found that finite-size effects for 5x5x5 
supercells are negligible for the quantities studied here 
(they are smaller than the statistical noise). 

In the path-integral formulation of statistical mechan- 
ics, the partition function is evaluated through a dis- 
cretization of the density matrix along cyclic paths, 
composed of a finite number TV-pr (Trotter number) of 
'imaginary-time' stepsii^ In the numerical simulations, 
this discretization gives rise to the appearance of Nt^ 
replicas for each quantum particle. In this way, the 



implementation of this method is based on an iso- 
morphism between the quantum system and a classi- 
cal one, obtained by replacing each quantum particle 
(here, atomic nucleus) by a cyclic chain of Nti clas- 
sical particles, connected by harmonic springs with a 
temperature-dependent constant. Details on this com- 
putational method can be found elsewhere 1^^*^^*^ 

To have a nearly constant precision for the simulation 
results at different temperatures, we considered a Trot- 
ter number that scales as the inverse temperature. At 
a given T, the actual value iVxr required to obtain con- 
vergence of the results depends on the Debye tempera- 
ture Qd (higher needs larger A^Tr)- For the simu- 
lations at zero pressure, we have taken N^^T — 250 K 
for solid Ar and N^^T = 200 K for the other rare-gas 
sohds {Qd ~ 90 K for Ar vs - 70 K for Ne, Kr, and 
Xe). Since vibrational frequencies (and the associated 
Debye temperature) increase as the applied pressure is 
raised, the Trotter number has to be correspondingly in- 
creased. Thus, for a given solid and an applied pressure 
we have taken iVxr values roughly proportional to the 
zero-point vibrational energy at the considered pressure. 
This means that N^^ is increased by a factor of about 
two for Ar, Kr and Xe (about three for Ne) when pres- 
sure rises from zero to 30 kbar. Thus, the computational 
time required to carry out PIMC simulations rises (a) as 
temperature is lowered (cc l/T), and (b) as pressure is 
raised [oc £'vib(0), zero-point vibrational energy]. For ex- 
ample, a PIMC simulation for solid Ar at 5 K and zero 
pressure {Nti = 50, N = 500) is equivalent in compu- 
tational time to a classical MC simulation for Nt^N = 
25000 atoms. This number increases by a factor of two 
at the same temperature and P = 30 kbar. 

Sampling of the configuration space has been carried 
out by the Metropolis method at temperatures between 
5 K and the triple-point temperature Ttp of the different 
solids, as well as at pressures up to 30 kbar. For given 
temperature and pressure, a typical run consisted of the 
generation of 2 x 10* quantum paths per atom for system 
equilibration, followed by 3 x 10^ paths per atom for the 
calculation of ensemble average properties. Ot her tec h- 
nical details are the same as those used in Refs. llSllOi 

The isothermal bulk modulus B can be obtained in 
the NPT ensemble from the mean-square fluctuations in 
the lattice parameter, cr^. In this ensemble, fluctuations 
in the volume V of the simulation cell are given by2i 



VkBT/B, and therefore 



B = 



(1) 



where L is the side length of the simulation cell in units 
of the lattice parameter (here, L = 5). 



B. Quasiharmonic approximation 

In the following section, results of PIMC simulations 
are compared with those derived from a QHA. This 
approximation is based on a renormalization of the 
phonon frequencies with volume, and for a given vol- 
ume the solid is assumed to be harmonic This vol- 
ume dependence of phonon frequencies is usually de- 
scribed by a mode-dependent Griineisen parameteriS 
7n(q) = — 91ncj„(q)/91ny, where a;„(q) are the fre- 
quencies of the nth mode in the crystal, and for small 
volume changes 7n(q) is assumed to be constant for each 
mode. However, for the QHA calculations presented 
here, we have not employed this description based on 
Griineisen parameters. Instead of this, we calculate di- 
rectly the actual (harmonic) vibrational frequencies for 
each crystal volume, by diagonalizing the corresponding 
dynamical matrix. 

For direct comparison with results of PIMC simula- 
tions, we have employed for the QHA the same supercell 
as for the simulations, i.e., a 5 x 5 x 5 supercell with 
periodic boundary conditions. This means that only the 
q = modes in the Brillouin zone of the supercell are in- 
cluded in the calculation, since modes with q 7^ violate 
the periodic boundary conditions. Then, the total num- 
ber of vibrational modes in the QHA is 1497, i.e., three 
times the number of rare-gas atoms in the supercell minus 
three translational degrees of freedom. The point group 
symmetry of the simulation cell imposes that many of 
these normal frequencies are degenerated. The number 
of normal modes that are not symmetry equivalent is 72 
for the supercell employed here. For each temperature, 
we calculated the free energy as a function of volume, 
with the corresponding phonon frequencies. The lattice 
parameter was changed in steps of 10"'^ A, and from the 
volume derivative of the free energy we derived the equi- 
librium volume as a function of pressure. 



III. RESULTS 
A. Energy 

Once defined an interatomic potential, the internal en- 
ergy of a solid, E{V, T), at given volume and temperature 
can be written as: 



EiV, T) = Eo + E,i{V) + E^,^{V, T) , 



(2) 



where Eq is the minimum potential energy for the (clas- 
sical) crystal at T = 0, Ec\{V) is the elastic energy, and 
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FIG. 1: Temperature dependence of the vibrational energy 
of solid neon. Circles and triangles correspond to kinetic and 
potential energy, respectively, as derived from PIMC simula- 
tions. Error bars of the simulation results are less than the 
symbol size. Dashed lines are guides to the eye. The solid 
line is the result of the QHA. Diamonds are kinetic-energy 
data for solid Ne, obtained by Timms et alw~ from neutron 
Compton scattering (the point at T = is an extrapolation 
given by these authors). 



_Bvib {y^ T) is the vibrational energy. Since we are work- 
ing in the isothermal-isobaric ensemble, it is understood 
that the volume is implicitly given by the applied pres- 
sure, i.e., V = V{P). For a given volume V , the classical 
energy at T = increases by an amount Ec\{V) with 
respect to the minimum energy Eq. This elastic energy 
Ec\ depends only on volume, but at finite temperatures 
and for the real (quantum) solids, it depends implicitly 
on T due to the temperature dependence of V (thermal 
expansion). The elastic energy E^iiV) represents a non- 
negligible part of the internal energy, even at zero pres- 
sure. For example, in Ne it is found to be 1.2 and 2.1 
meV per atom at 5 and 24 K, respectively. These values 
are smaller for the other rare gases, as shown in Table I. 

The vibrational energy, _Evib(^, ?^), depends explicitly 
on both, V and T, and can be obtained by subtract- 
ing the elastic energy from the internal energy. Values 
of i?vib derived from our PIMC simulations for T ^ 
and P = are given in Table I. Path-integral Monte 
Carlo simulations allow us to obtain separately the ki- 
netic, and potential energy, iJp, associated to the 
lattice vibrationsi^ Both energies are shown in Fig. ^for 
solid Ne as a function of temperature at P = 0. Circles 
and triangles correspond to the vibrational kinetic and 
potential energy, respectively. Our results for the kinetic 
energy are close to those derived earlier from PIMC simu- 



lations with Lennard-Jonesi^iSii^ and Azia^S interatomic 
potentials. For comparison, we present also in Fig. Rval- 
ues of the kinetic energy of Ne atoms, derived by Timms 
et al/^^ from neutron Compton scattering in solid neon 
(black diamonds). According to the results of our PIMC 
simulations, Ek is larger than Ep by about 20%. The 
QHA predicts potential (and kinetic) energy values (solid 
line) which are close to the vibrational potential energy 
derived from our PIMC simulations. Something similar 
happens for the other rare-gas solids, with E^^ > Ep for 
all temperatures and pressures studied here. The differ- 
ence Ek — Ep decreases for increasing atomic mass, and 
at T = 5 K and zero pressure, we find Ek — Ep = 0.67 
and 0.059 meV/atom for Ne and Xe, respectively. These 
energy differences increase slowly as pressure rises, and 
take values of 0.75 and 0.060 meV/atom for Ne and Xe 
at 30 kbar. 

A quantitative estimation of the overall anharmonic- 
ity of the atom vibrations is given by the parameter—: 
i — 2{Ek — Ep)/{Ek+Ep), which should be zero for a har- 
monic solid at any temperature, as follows from the virial 
theorem. For rare-gas solids, it was shown earlier— that 
^ increases as temperature rises, as expected for larger 
anharmonicity. In Fig. [21 we show the pressure depen- 
dence of the parameter ^ for different rare-gas solids at 
T = 5 K, as derived from PIMC simulations. One ob- 
serves that ^ decreases as pressure is raised. The relative 
change in this parameter is largest for Ne, for which it 
decreases by a factor of about 3.5. For Xe, it changes by 
a factor « 2. 

The elastic energy Ec\ increases fast as pressure rises. 
In Fig. 121 we display the pressure dependence of E^i for 
(a) Ne and (b) Ar, obtained from PIMC (open squares) 
and QHA (solid lines). For comparison we also present 
results for the vibrational energy (kinetic plus potential) 
obtained by both procedures. As shown above, Ek for Ne 
at zero pressure is about 20% larger than Ep, and thus 
the whole i?vib derived from PIMC is 10% larger than 
that found in the QHA {Ep is nearly the same for both 
methods). This relative difference decreases as pressure 
increases, and is almost inobservable on the scale of Fig. 
O Note that close to P = 0, E^ decreases as P rises, 
reaches a minimum and then it increases continuously 
for larger applied pressures. This is due to the fact that 
at P = the crystal is expanded with respect to the 
volume giving the minimum potential energy. With an 
applied pressure of about 2 kbar, solid neon reaches the 
volume corresponding to the classical minimum (giving 
£-01 = 0), and for larger pressures the volume is further 
reduced, with an increase in Ed- Something similar hap- 
pens for Ar and the other rare-gas solids, although it is 
less appreciable on the scale of Fig. |SIb). For solid Ne 
we find that the vibrational energy is larger than E^i at 
pressures P < 23 kbar. However, E^i rises faster with 
pressure, and becomes the dominant part of the inter- 
nal energy for larger pressures. This behavior is quali- 
tatively similar for the other rare-gas solids, as shown in 
Fig. Olb) for Ar. The main difference is that the elas- 
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FIG. 2; Anharmonicity parameter ^ for rare-gas crystals as 
a function of applied hydrostatic pressure. Symbols indicate 
results of PIMC simulations at T = 5 K. From top to bottom: 
Ne, Ar, and Xe. Error bars are less than the symbol size. 
Lines are guides to the eye. Results for Kr lie between those 
for Ar and Xe, and are not shown for clarity of the figure. 



tic energy for Ar increases with pressure faster than for 
Ne, and becomes larger than i?vib for P > 12 kbar. In 
this context, the main point concerning large pressures 
is that dEc\/dP > dEvih/dP. This means that the ra- 
tio Eci/Evih grows for increasing pressure, and eventually 
the vibrational energy becomes a small correction to the 
internal energy. 

Note that in Fig. Owe have included some points at 
negative pressure (that is, solids under tension). This re- 
gion with P < was studied earlier for rare-gas solids by 
PIMC simulations^ It was found, in particular, that at 
T = 5 K solid Ne and Ar are metastable until reaching the 
corresponding spinodal pressures of -0.9 and -2.5 kbar, 
where they become mechanically unstable. Here we only 
comment that energy results obtained from the QHA fol- 
low those yielded by PIMC simulations also in this region 
of negative pressures. 



B. Heat capacity 

We have calculated the heat capacity Cp of rare-gas 
solids for several pressures as a numerical derivative of 
the enthalpy, H — E + PV, with respect to the temper- 
ature. In Fig. ^we present results derived from PIMC 
simulations (symbols) for solid argon at P = and 15 
kbar. Results of the quasiharmonic approximation are 
shown as dashed lines. For comparison, we also present 
experimental results obtained by Flubacher et al^ for 
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FIG. 3; Vibrational and elastic energy of solid Ne (a) and 
Ar (b) at T = 5 K, as a function of pressure. Symbols show 
results of PIMC simulations: squares, vibrational energy; cir- 
cles, elastic energy. Solid lines are results of the quasihar- 
monic approximation. Error bars of the simulation results 
are less than the symbol size. 



argon at atmospheric pressure. Results of PIMC simu- 
lations at P = agree well with experimental data for 
both Ne and Ar (data for Ne were given elsewherai^ and 
are not shown here), except close to Ttp, where the sim- 
ulation results are lower than the experimental ones in 
both cases. However, for T < Ttp there can be a sys- 
tematic error in the experimental results ('^ 10% at Ttp), 
due to partial vaporizing of the solid4i Our simulations 
were carried out for perfect crystals (without vacancies), 
and thus we cannot conclude at this point if the observed 
difference is caused by a failure of the Lennard-Jones po- 
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FIG. 4: Heat capacity, Cp, of solid argon as a function of 
temperature at P = and 15 kbar. Symbols: PIMC sim- 
ulations; dashed lines: quasiharmonic approximation; solid 
line: experimental data obtained by Flubacher et a/.''° at at- 
mospheric pressure. Error bars of the simulation results for 
P = are on the order of the symbol size. 



tential employed here, or by the experimental uncertainty 
near Ttp. 

For zero pressure, the QHA gives results close to those 
of PIMC at low temperatures (T < 5 K for Ne and 
T < 20 K for Ar), but clearly overestimates the heat 
capacity at higher temperatures. The relative error of 
the QHA (compared with PIMC) increases as tempera- 
ture is raised, as expected for an increase in (intrinsic) 
anharmonicity due to thermal effects. In the case of Ne, 
this approximation breaks down at T = 22.6 K, due to a 
mechanical instability (the solid becomes unstable in the 
QHA at a temperature lower than the actual Ttp). For 
Ar this approach predicts a breakdown of the solid at 
86.2 K. These results are in line with earlier observations 
on the validity temperature range of the QHA for these 



For increasing pressure, the error of the QHA de- 
creases. This is in part due to the renormalization of 
vibrational frequencies, which increase under an applied 
pressure, and therefore shift the increase in Cp to higher 
temperatures. Moreover, the relative contribution of the 
vibrational energy to the internal energy is reduced as 
pressure rises, and in consequence the intrinsic anhar- 
monicity of the vibrational modes (not captured by the 
QHA) becomes less relevant for the heat capacity. 



C. Lattice parameter 

At T = 0, the difference Aa(0) = a(0) - aci(O) be- 
tween tfie actuai iattice parameter a(0) and that corre- 
sponding to the minimum potential energy of the (clas- 
sical) crystal, aci(O), decreases as the atomic mass rises 
and quantum effects become less relevant»ii^ From our 
PIMC simulations we found at P = that Ao(0) ranges 
from 0.174 A for Ne to 0.025 A for Xe. Calculated val- 
ues for a(0) at zero pressure are given in Table I. The 
temperature dependence of the lattice parameter derived 
from this kind of PIMC simulations agrees with exper- 
imental data for rare-gas solidsji2ii2i Such an agreement 
is also found at low temperatures in the pressure range 
considered here (P < 30 kbar). In this pressure range, 
the QHA predicts lattice parameters which follow closely 
those given by the simulations (see below the discussion 
on the bulk modulus). 

For large pressures, it is known that the description 
of rare-gas solids with effective interatomic potentials, 
requires the consideration of three-body termsi^iS^i^ to 
reproduce the actual equation of state P-V. We have 
checked that the Lennard-Jones potential considered here 
predicts P-V isotherms for solid Ar close to the exper- 
imental ones2^^ up to pressures on the order of 50 kbar. 
For larger pressures, this pair potential yields volumes 
larger than the real ones. 

Even though thermal effects on the crystal volume 
change in magnitude for different rare-gas solids at low 
temperatures (T < 5 K), these differences are less im- 
portant at higher temperatures. If one takes as a refer- 
ence the classical lattice parameter aci(O), the difference 
Aa = a — aci(O) at temperatures close to Ttp of each solid 
amounts to « 0.25 A, as derived from our PIMC simu- 
lations. This difference Aa is plotted in Fig. [S] for the 
different solids as a function of pressure. For each solid 
we present PIMC results along an isotherm: Ne, T = 20 
K; Ar, 80 K; Kr, 110 K; Xe, 160 K. These results for the 
different rare-gas solids follow roughly the same pressure 
dependence. Note that at P = 0, Aa is positive, as a 
consequence of zero-point and thermal lattice expansion. 
Under applied pressure Aa decreases, reaching Aa = 
for a pressure P w 2 kbar, and is negative for larger P. 
Lines in this figure show results of the QHA, which follow 
closely those of PIMC. The main difference between both 
sets of results appears for Xe at pressures higher than 10 
kbar (diamonds and dashed-dotted lines). We conclude 
that the difference Aa at temperatures near Ttp follows 
a pressure dependence similar for all rare-gas solids con- 
sidered here. 



D. Bulk modulus 

The isothermal bulk modulus B of rare-gas solids has 
been calculated from our PIMC simulations by using 
Eq. The temperature dependence of B at zero pres- 
sure was studied earlieriiSiiS and will not be presented 
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FIG. 5: Pressure dependence of the increment in lattice pa- 
rameter, Aa — a~ aci(O), with respect to the ideal (classical) 
crystal with minimum potential energy. Symbols and lines 
indicate results of PIMC simulations and QHA along a given 
isotherm for each solid: Ne at 20 K (squares, solid line), Ar 
at 80 K (circles, dotted line), Kr at 110 K (triangles, dashed 
line), and Xe at 160 K (diamonds, dashed-dotted line). 



here. We only note that results obtained from this kind 
of simulations for Ar, Kr, and Xe showed good agree- 
ment with experimental data up to temperatures close 
to TtpiiSiiS Results found for the bulk modulus of neon 
are lower than the experimental data at T < 18 K, and at 
higher T the experimental B decreases faster than that 
derived from the simulations. This seems to be a general 
problem of this kind of calculations with Lennard-Jones 
and Tang-Toennis-type interatomic potentialsiSSii^ 

We now turn into the pressure dependence of the bulk 
modulus for these solids. Results of our PIMC simula- 
tions for Ne, Ar, and Xe at T = 5 K are plotted in Fig. 
El (symbols). The bulk modulus of Kr (not shown for 
clarity) lies between those of Ar and Xe. Dashed lines 
are results of the QHA, and solid lines represent exper- 
imental data obtained by Anderson et alM- for Ne, and 
Anderson and Swenson'*^ for Ar and Xe. The QHA pre- 
dicts bulk moduli in good agreement with those derived 
from PIMC simulations, indicating that this approxima- 
tion is rather accurate for predicting the P-V equation 
of state, as well as the derivative dP/dV, which gives the 
bulk modulus. 

Our calculated results for Ne agree well with experi- 
mental data in the pressure region under consideration. 
For the other rare-gas solids (including Kr, not shown in 
Fig. our results are larger than the experimental data 
for P > 10 kbar. At 20 kbar, our method overestimates 
the bulk modulus of Ar and Xe by about 5%. This is a 
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Pressure (kbar) 

FIG. 6: Isothermal bulk modulus as a function of pressure for 
rare-gas solids. Symbols: results of PIMC simulations at T = 
5 K, with error bars on the order of the symbol size. Dashed 
lines: quasiharmonic approximation at 5 K. Solid lines: ex- 
perimental data by Anderson et al^^ for Ne, and Anderson 
and Swenson— for Ar and Xe, at T = 4.2 K. 



consequence of the interatomic pair potential employed 
in our calculations. As indicated above, the limitations 
of pair potentials for describing this kind of solids show 
up as pressure rises, and eventually three-body terms are 
necessary at high pressures li^i^ 



IV. DISCUSSION 

Path-integral Monte Carlo simulations give in principle 
the "exact" solution to the considered quantum problem, 
with an accuracy depending on the considered Trotter 
number and the statistical error associated to the MC 
sampling. Thus, the Lennard-Jones potential employed 
here gives a good description of structural and thermo- 
dynamic properties of rare-gas solids at zero pressure be- 
tween T = and the triple-point temperature. The equa- 
tion of state P — V is well described by this potential in 
the pressure range considered here (P < 30 kbar). The 
bulk modulus of Ne is well reproduced, and for heavier 
rare gases PIMC simulations give B values larger than 
experiment at P > 10 kbar. 

In all solids considered here, and for different pressures 
and temperatures, we have found Ek > Ep. The overall 
anharmonicity has been measured by the parameter ^, 
which for given T and P decreases as the atomic mass is 
raised. Although this parameter is by no means a unique 
and absolute measure of the whole anharmonicity, for a 
family of similar materials it is useful to give us a quan- 



titative estimation of anharmonic effects as a function 
of atomic mass, temperature, and pressure. In the limit 
P = and T ^ 0, ^ changes from 0.21 for Ne to 0.02 for 
Xe. This indicates a large anharmonicity of the lattice 
vibrations in Ne, even at T = 0. For comparison, we note 
that covalent solids at low temperatures show much lower 
values of ^. Thus, for diamond, silicon, and germanium 
one finds for T/Qd ^ 1, differences between Ep and 
Ek smaller than 1%, and aXT ^ Qo they are less than 
3%4ii2ii£ Another point of interest is that for such cova- 
lent materials the vibrational Ep was found to be larger 
than Ek, just the opposite to the trend found for rare- 
gas solids. Then, the fact that Ep < Ek, obtained for the 
rare-gas solids studied here, is not general in solids, and 
can be due to the particular nature of the interatomic 
interactions present in rare-gas (Lennard-Jones) solids. 

A qualitative understanding of the sizeable increase in 
kinetic energy of the rare-gas atoms with respect to the 
value expected in a QHA can be obtained by analyzing 
the changes of kinetic and potential energy by standard 
time-independent perturbation methods. With this pur- 
pose, we consider a one-coordinate perturbed harmonic 
oscillator with Hamiltonian H = /2m + ^muPx^ -I- 
Ax^ -t- Bx"^. The x^ term does not introduce corrections 
to the zero-point energy in first order, '^^ and the second- 
order correction is only due to a change in kinetic energy. 
The x'^ term gives a first-order correction^ that again 
is only caused by a kinetic- energy change. Thus, the 
leading corrections to the zero-point energy due to both 
perturbing terms originate from changes in the kinetic en- 
ergy. This is in line with the results presented in Fig. ^ 
which show that the potential energy derived from PIMC 
is close to that yielded by the QHA, and the kinetic en- 
ergy is far from the QHA result. This one-coordinate 
approach can give us only a very qualitative interpre- 
tation of the changes in kinetic and potential energy of 
the considered crystals with respect to a quasiharmonic 
approximation. In fact, the whole problem is a many- 
particle one, in which the lattice vibrations cannot be 
considered as non-interacting entities when anharmonici- 
ties are present. The coupling between vibrational modes 
obtained in a harmonic or quasiharmonic approximation 
is expected to increase as the crystal density increases 
(pressure rises), thus making such approximations less 
reliable in the description of the vibrational problem. 

For what concerns structural and thermodynamic 
properties of rare-gas solids, the QHA becomes more pre- 
cise as pressure increases. This is related to the pressure 
dependence of the vibrational and elastic energies, and 
is due to two reasons. First, for high pressures the vi- 
brational energy becomes a small part of the internal en- 
ergy, and hence the influence of lattice vibrations (and of 
their anharmonicity) on thermodynamic properties will 
comparatively decrease. The same happens for the free 
energy F at finite temperatures, since its vibrational part 
Pvib becomes small as compared with the whole F . Thus, 
for Pvib ^ F — Eq, a QHA describes well the pressure de- 
pendence of the crystal volume, since it is basically given 
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by the elastic energy of the solid. Second, the intrin- 
sic anharmonicity of the lattice vibrations (as measured 
by the anharmonicity parameter ^) decreases as pressure 
rises (see Fig. and eventually becomes negligible for 
large P (^ ^ 0). Both arguments go in the same direc- 
tion of making more accurate the QHA. 

Nevertheless, the improved accuracy of the QHA as 
pressure rises is not a particular merit of this approach, 
since the internal energy becomes dominated by the elas- 
tic energy and the actual description of the vibrational 
modes is not very relevant for thermodynamic proper- 
ties. It has been also arguedS that structural properties 
of solids under pressure can be described rather accu- 
rately by a classical model for the lattice vibrations. The 
origin of this is similar to that described above for the 
decreasing effect of anharmonicity as P rises, since in this 
respect the actual description of the lattice vibrations by 
a classical or a quantum model becomes unimportant for 
solids under large pressures. This is not necessarily the 
case for spectroscopic properties of the solids under con- 
sideration, because vibrational frequencies predicted by 
a QHA or by a classical model are not guaranteed to 
describe correctly the actual ones for high P. 

In summary, we have analyzed the influence of a hydro- 
static pressure on anharmonic effects of rare-gas solids. 
Our results indicate that the validity of the QHA to de- 
scribe structural and thermodynamic properties of these 
solids increases as pressure is raised. This is mainly a 
consequence of the relative importance of elastic and vi- 
brational energy, since the latter becomes increasingly 



irrelevant as pressure rises. Therefore, the precision re- 
quirements for a description of the (anharmonic) vibra- 
tional modes is reduced with increasing pressure. In the 
limit of large pressures, even a classical description of 
these modes can be sufficiently precise to predict several 
properties of these solids at low temperatures. This is, of 
course, not the case of vibrational properties, which may 
require the full quantum treatment with the considera- 
tion of zero-point anharmonic effects at high pressures. 

We finally note that the extension of the method em- 
ployed here to study solids under larger pressures, such as 
those presently reached in experimental studies [P ^ 100 
GPa) , is hampered by the requirement of enlarging enor- 
mously the Trotter number (and consequently the CPU 
time) in PIMC simulations, which has to be increased 
as the vibrational energy rises. This problem is partic- 
ularly important to study vibrational properties, since 
for very high pressures, thermodynamic properties can 
be well described by neglecting the quantum nature of 
the atomic nuclei, as usually done in electronic-structure 
calculations. 
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